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1.  Introduction 

Recent  and  projected  developments  in  supercomputers,  numerical  grid  gen¬ 
eration  techniques  and  computational  algorithms  for  the  compressible  Euler  and 
Navier-Stokes  equations  portend  a  major  revolution  in  the  manner,  pace,  and  cost 
of  design  and  the  resulting  performance  of  aerodynamic  systems.  To  realize  these 
potential  benefits,  certain  closing  developments  in  computational  technique  must 
be  made  in  order  to  effect  a  highly  accurate,  reliable,  efficient  and  productive  sim¬ 
ulation  environment  for  aerodynamic  design  analysis. 

A  primary  need  of  the  developments  is  to  achieve  the  capability  for  a  user  to 
easily,  rapidly  and  accurately  perform  flowfield  calculations  among  problems  of  dis¬ 
parate  and  realistically  complex  geometries.  The  natural  approach  to  realizing  this 
objective  with  comparatively  straightforward  extensions  of  existing  finite  difference 
computational  technology  is  through  the  use  of  systems  of  quadrilateral  patched 
meshes. 

Such  systems  can  be  either/both  composite  (joined)  or  overset  (disjoint).  In 
the  former  case  adjacent  patches  share  a  common  boundary,  or  at  least  parallel 
boundaries  in  the  case  of  mesh  patch  overlap  for  purposes  of  applying  numerical 
boundary  conditions.  With  composite  meshes,  patch  boundaries  are  piecewise  fitted 
to  segments  of  physical  or  computational  boundaries  or  embedded  flow  structures. 
As  shown  by  Lombard,  et  al1,  composite  mesh  systems,  that  may  have  numerically 
useful  properties  of  geometric  continuity  across  patch  boundaries,  admit  topolog¬ 
ically  singular  global  meshes  that  have  the  capability  to  connect  computational 
regions  of  great  (really  any)  geometric  complexity.  However,  situations  exist  where 
multiple  mesh  topologies,  each  naturally  related  to  some  different  piece  of  geome¬ 
try  or  flow  structure,  offer  greater  flexibility  and  accuracy  than  composite  meshing 
alone.  Examples  involve  multiple  bodies  that  may  have  relative  motion  or  weak 


shocks  that  bear  little  geometric  relationship  to  boundaries  of  the  flow  domain.  In 
such  cases  systems  of  disparately  oriented  and  at  least  partially  overset  grids  as 
proposed  by  Berger  and  Oliger2  allow  arbitrarily  high  resolution  of  all  features  of 
the  flowfield. 

To  make  efficient,  productive  use  of  patched  meshing  strategies  requires  a  body 
of  new  computational  tools  and  methodology  that  are  the  objectives  of  the  present 
research.  The  needed  factors  are:  (1)  a  simple  procedure  for  generating  patched 
computational  meshes  with  freedom  of  point  and  gradient  specification  on  all  patch 
boundaries;  (2)  improved  upwind  algorithm/numericai  boundary  condition  proce¬ 
dures  for  semi-autonomous  implicit  but  unconditionally  stable  conservative  coupling 
of  solutions  on  a  system  of  multiple  patched  meshes;  and  (3)  computer  graphics, 
particularly  a  simple  algorithm  for  constructing  contour  plots  on  systems  of  over¬ 
set  patched  meshes.  The  glue  that  is  to  tie  these  tools  together  in  a  simulation 
requiring  minimum  human  intervention  to  adapt  to  new  configurations  is  a  flexible 
parameter  controlled  multiple  mesh  data  structure.  An  important  objective  of  the 
program  is  to  test  the  evolving  techniques  in  appropriate  problems. 

2.  Research  Accomplishments 

This  report  concludes  contract  F49620-83-C-0084  and  summarizes  the  research 
accomplishments  of  the  first  two  years  of  an  anticipated  three  year  program.  In  the 
first  year  the  emphasis  was  placed  on  the  most  crucial  and  challenging  objectives 
-  patched  grid  generation  and  robust  upwind  algorithm/boundary  procedures  for 
rapid  relaxation  on  multiple  meshes.  The  first  years  effort3  sufficed  to  create  some 
needed  tools  of  algebraic  grid  generation  and  to  establish  an  operationally  explicit 
but  unconditionally  stable  upwind  algorithm/numerical  boundary  condition  proce¬ 
dure  for  systems  of  patched  meshes.  The  techniques  were  tested  in  comparatively 
simple  problems. 

In  the  second  years  effort  the  techniques  were  tested  against  challenging  inter¬ 
nal  aerodynamics  problems  and  improved  and  extended  from  two  to  three  space 
dimensions.  A  simple,  effective  algorithm  for  contour  plotting  on  domains  covered 
by  multiple  patched  meshes  was  created  to  exhibit  results. 

Algebraic  Grid  Generation 

The  concept  of  patched  meshing  in  which  complex  domains  are  broken  up  into 
many  geometrically  regular  and  topologically  rectangular  subdomains  leads  natu- 


rally  to  the  use  of  efficient  algebraic  techniques  for  the  construction  of  the  individual 
mesh  patches.  To  obtain  the  desired  smoothness  properties  over  the  global  mesh 
in  the  vicinity  of  patch  boundaries,  a  technique  that  permits  specification  of  point 
distribution  and  gradient  on  all  boundaries  was  devised.  The  technique  -  termed 
generalized  transfinite  interpolation4  -  makes  use  of  a  parameterized  general  cubic 
polynomial  for  the  coordinate  curves.  Regularity  of  the  mesh  is  obtained  by  em¬ 
ploying  continuous  distributions  of  the  parameters  of  the  curves  within  judiciously 
chosen  bounds  based  on  analysis.  Stretching  functions  such  as  that  of  Vinokur5  are 
used  to  distribute  points  and  blending  functions  are  used  to  distribute  parameters 
of  the  curves  between  lateral  boundaries. 

A  novel  feature  of  the  technique  is  the  introduction  of  the  corner  singularity 
from  analysis  to  govern  distribution  of  points  and  parameters  in  the  vicinity  of 
boundary  slope  singularities.  At  such  points,  the  method  thus  obtains  the  desired 
properties  of  mesh  smoothness  to  the  interior.  A  global  mesh  solution  obtained  with 
the  method  for  a  backward  step  problem  is  shown  in  Figure  1.  Here  the  solution 
was  generated  in  two  patches  one  containing  the  exterior  corner  and  the  other  the 
interior  corner.  The  solution  was  matched  analytically  at  the  patch  interface. 

Early  in  the  second  year,  in  the  process  of  attempting  to  apply  the  generalized 
transfinite  interpolation  technique  in  a  variety  of  2-D  problems  it  became  evident 
the  method  was  too  sensitive  to  parameter  selection  among  too  many  options,  was 
confusing  and  ultimately  required  too  much  artistry  to  meet  the  objectives  of  sim¬ 
plicity  and  user  friendliness  set  for  the  products  of  the  research.  Further,  the  lack 
of  an  analytical  solution  to  corner  problems  blocked  the  straightforward  extention 
of  the  technique  to  3-D. 

With  some  reflection  it  became  evident  to  us  that  the  difficulty  lay  in  trying 
to  accomplish  too  much  in  a  single  step  process.  Rather,  borrowing  the  tools  of 
the  algebraic  technique  and  redefining  the  process  in  multiple  steps  with  interactive 
computer  graphics  sets,  we  could  define  a  straightforward  procedure  to  meet  the 
desired  ends. 

The  approach  that  has  been  implemented  is  in  the  realm  of  two  boundary 
methods  in  that  one  pair  of  opposite  sides  of  a  patch  is  regarded  as  prescribed  and 
often  includes  a  portion  of  a  physical  boundary.  The  other  pair  of  sides  is  formed 
of  the  left  and  right  limiting  members  of  the  family  of  generalized  cubic  coordinate 


curves  joining  the  initially  given  two  boundaries.  In  either  2-D  or  3-D  the  general 
cubic  coordinate  curve  has  the  simple  form 


t  =  £1  +  (r?  -  rjJ/(u)  +  aig{u)  +  ^(u) 


(1) 


where 

f(u )  =  u2(3  -  2u) 
g(u)  =  u(l  -  u )2 
h(u)  —  u2{u  -  1) 

Equation  (l)  is  a  hermite  interpolation  of  value  (r)  and  gradient  (er)  on  the  two 
boundaries  and  is  parameterized  in  terms  of  u  which  varies  from  zero  to  unity.  The 
scalings  of  q_ j_  and  02  influence  the  shape  (curvature)  of  the  curve  between  any  pair 
of  end  points.  The  specification  of  a  discrete  set  of  u  values  using  a  generalized 
distribution  function  such  as  that  of  reference  4  defines  the  nodal  intersections  with 
the  other  family  of  coordinate  curves. 

In  our  implementation  the  left  and  right  limiting  (lateral  bounding)  coordinate 
curves  are  developed  interactively  on  a  graphics  terminal  or  workstation  to  have 
the  desired  configurations.  The  parameters  of  these  lateral  bounding  curves  are 
then  blended  with  polynomial  weighting  functions  to  describe  the  general  cubic 
coordinate  curve  over  the  intervening  also  discretized  interval. 

The  lateral  patch  boundaries  are  essentially  control  devices  that  specify  shape 
and  distribution  to  surrounding  regions.  As  such  they  are  placed  where  needed  - 
at  breaks  in  body  surface  geometry  and  as  terminators  or  transition  guides  from 
regions  of  strong  shape  variation  to  regions  of  very  regular  mesh.  In  fact  the  mesh 
generation  problem,  particularly  for  geometries  with  any  substantial  complexity,  is 
a  problem  of  multiple  length  scales.  The  purpose  of  multiple  patching  is  to  isolate 
regions  of  comparable  scales  and  on  which  subdomains  the  solution  is  comparatively 
regular  and  can  be  conveniently  fit  by  simple  functions. 

Once  a  primary  grid  is  generated  by  the  technique  described  above  it  can  be 
interactively  improved  by  modifying  parameter  blending  and  point  distribution  in¬ 
cluding  point  redistribution  along  the  alternate  family  of  coordinate  lines  implicitly 
defined  by  the  nodes  on  the  cubic  coordinate  curves.  The  latter  operation  is  in 
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the  spirit,  if  not  the  detailed  implementation  of  a  two  step  generalized  transfinite 
interpolation. 

Another  secondary  operation  that  we  employ  is  the  modification  of  coordinate 
lines  in  the  vicinity  of  a  boundary  to  smoothly  enforce  local  normality.  The  latter 
operation  like  all  the  procedures  has  been  programmed  as  a  convenient  tool  requiring 
minimal  input  to  apply  at  a  boundary.  Finally,  a  parameterized  tension  spline  that 
provides  an  analytical  description  of  a  curve  amongst  discrete  data  is  a  tool  that 
has  proved  useful  in  the  latter  operation,  in  effecting  redistribution  of  points  along 
coordinate  curves  of  either  family  and  for  fitting  numerically  specified  boundary 
data.  In  Figures  2  and  3  are  shown  respectively  typical  examples  of  2-D  and  3-D 
mesh  planes  generated  with  the  simplified  technique. 

Upwind  Implicit  Relaxation  Algorithm/Boundary  Procedures 

Under  the  contract  we  have  devised  a  new  single  level  operationally  explicit 
but  effectively  implicit  algorithm  for  gasdynamics.  The  algorithm  is  particularly 
appropriate  for  multiple  patch  mesh  systems  because  each  solution  sweep  operation 
on  any  patch  is  decoupled  from  any  other.  Thus  the  method  is  not  only  very  storage 
efficient  and  simple  to  program  including  the  coupling  at  patch  boundaries  but,  also, 
can  make  excellent  use  of  parallel  computing  in  several  straightforward  ways. 

Previously  the  Beam- Warming  factored  implicit  algorithm6  with  the  Baldwin- 
Lomax  thin  layer  viscous  approximation7  has  provided  the  basis  for  two  similar 
space  marching  (PNS)  procedures8,9  for  the  compressible  Navier-Stokes  equations. 
These  PNS  methods  which  are  highly  efficient  -  requiring  half  the  data  storage 
and  a  small  fraction  of  the  computer  time  of  two  level  time  dependent  methods  - 
have  proven  effective  for  flows10,11  with  favorable  streamwise  pressure  gradient  or 
with  relatively  small  adverse  pressure  gradients.  However,  in  the  presence  of  strong 
adverse  pressure  gradient  such  as  occurs  in  a  wing  or  fin  root  regions  the  contem¬ 
porary  PNS  methods  suffer  numerical  stability  problems  and  may  infer  streamwise 
separation  even  where  separation  doesn’t  occur12.  In  such  unseparated  (perhaps 
weakly  separated)  regions,  numerical  stability  has  been  maintained  at  the  price  of 
employing  large  amounts  of  artificial  viscosity  with  a  resulting  loss  in  predictive 
accuracy  and  knowledge  of  the  actual  state  of  the  flow.  Where  strong  streamwise 
separation  occurs  the  methods  are  unstable  and  cannot  proceed.  Particularly  for 
the  increasingly  relevant  laminar  flow  situation  that  will  be  encountered  at  very  high 
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altitude  by  aerodynamic  systems  such  as  orbital  transfer  vehicles  (OTV’s)  and  space 
shuttle,  streamwise  separation  becomes  a  likely  occurrence13  in  compression  corners 
associated  with  canopies,  pods,  flared  bodies,  wing  or  fin  roots  and  deflected  con¬ 
trol  surfaces.  Thus  a  more  general  technique  is  needed  that  is  inherently  stable  for 
all  types  of  upstream  influence.  At  a  minimum  the  mixed  elliptic-hyperbolic  prob¬ 
lem  requires  global  iteration,  preferably  with  type  dependent  differencing.  More 
background  on  this  problem  area  is  given  in  the  first  annual  report3. 

New  Universal  Single  Level  Scheme  CSCM-S 

The  CSCM  flux  difference  eigenvector  split  upwind  implicit  method  14'15'1G  for 
the  inviscid  terms  of  the  compressible  Navier-Stokes  equations  provides  the  natural 
basis  for  an  unconditionally  stable  space  marching  technique  through  regions  of 
subsonic  and  streamwise  separated  flow.  In  such  regions  the  split  method  can  be 
likened  to  stable  marching  of  each  scalar  characteristic  wave  system  in  the  direction 
of  its  associated  eigenvalue  (simple  wave  velocity).  In  supersonic  flow,  where  all 
eigenvalues  have  the  same  sign,  the  method  automatically  becomes  equivalent  to  the 
referenced  PNS  techniques  based  on  the  Beam- Warming  factored  implicit  method 
with  the  Baldwin-Lomax  thin  layer  viscous  approximation. 

Compared  to  contemporary  central  difference  methods,  the  CSCM  character¬ 
istics  based  upwind  difference  approximation  with  its  inherent  numerical  stability 
leads  to  greatly  reduced  oscillation  and  greater  accuracy  in  the  presence  of  captured 
discontinuities  such  as  shocks,  contacts  and  physical  or  computational  boundaries. 
The  correct  mathematical  domains  of  dependence  that  correspond  with  physical 
directions  of  wave  propagation  are  coupled  with  well  posed  characteristic  boundary 
approximations13  naturally  consistent  with  the  interior  point  scheme.  The  result 
is  faster  sorting  out  of  transient  disturbances  and  substantially  more  rapid  conver¬ 
gence  to  the  steady  state.  The  splitting  and  the  associated  time  dependent  implicit 
method  have  been  described  in  detail  in  references  (14)  and  (16)  for  quasi  1-D  and 
2-D  planar  or  axisymmetric  flow. 

In  the  following,  we  will  sketch  the  differences  between  the  time  dependent 
method  and  the  new  space  marching  technique  which  we  designate  CSCM-S.  The 
discussion  will  begin  with  the  quasi  1-D  inviscid  formulation,  present  some  results 
elucidating  the  properties  and  performance  of  the  method,  then  give  additional  de¬ 
tails  entering  into  multidimensional  inviscid  and  thin  layer  viscous  procedures  and 
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present  some  2-D  solutions  obtained  with  the  new  single  level  scheme  in  problems 
solved  previously16  with  the  time  dependent  method.  Lastly  we  sketch  a  3-D  im¬ 
plicit  method  of  planes  algorithm  and  give  some  results  for  an  axisymmetric  nozzle 
flow  over  a  backward  facing  step. 

Quasi  1-D  Formulation 

The  general  jth  interior  point  difference  equations  for  the  time  dependent 
CSCM  upwind  implicit  method  is  written 

{I  +  A+V  +  A-A)8q,- = -A+Aq)  U  -A~Aq)n  (2) 

3  ~  1  3 

where  V  and  A  are  backward  and  forward  spatial  difference  operators.  In  the 
notation  the  interval  averaged  matrices  between  node  points  j  and  j  + 1  are  indexed 
j.  The  right  hand  side  of  equation  (2)  is  written  for  the  first  order  method.  Higher 
order  methods  in  space  are  given  with  results  in  references  14  and  16. 

Central  to  its  accurate  shock  capturing  capability,  the  CSCM  conservative  flux 
difference  splitting  has  the  “property  U”  put  forth  by  Roe17 

U+  +  A~)Aq)3-  =  AF)y  =  Fj+1  -  Fj  (3) 

Here  q  is  the  conservative  dependent  variable  vector  and  F  is  the  associated  flux 
vector.  The  matrices  A+  and  A~  are  the  splittings  of  the  CSCM  interval  averaged 
Jacobian  matrix  according  to  the  signs  of  the  averaged  eigenvalues.  Thus  in  the 
equation  for  the  yth  grid  point,  A+ A^g)j_i  represents  stable  characteristic  spatial 
differencing  backward  for  positive  eigenvalue  contributions  and  A~A^q):i-,  forward 
for  negative  ones. 

With  Sq  =  qn+1  —  qn,  equation  (2)  defines  a  two  level  linearized  coupled  block 
matrix  implicit  scheme  that  can  be  solved  by  a  block  tridiagonal  procedure.  In 
reference  (16)  a  new  (DDADI)  approximately  factored  alternating  sweep  bidiagonal 
solution  procedure  for  equation  (2)  is  presented  that  is  shown  to  be  very  robust 
and  is  operationally  explicit,  i.e.  requires  only  a  decoupled  sequence  of  local  block 
matrix  inversions  rather  than  the  solution  of  the  coupled  set.  For  the  forward  sweep 
the  bidiagonal  solution  procedure  can  be  written 


For  the  linear  problem,  i.e.  constant  coefficient  case  of  stability  analysis,  equation 
(4)  is  equivalent  to  the  single  level  space  marching  procedure 


(/  +  A+  -  A-)6q*s  =  A+q  ,  J  -  A  +  qn  -  A~ Aq)  ”  (5) 

3  3 

Nonlinearity  enters  in  the  single  level  space  marching  form  (5)  in  that  at  each  step  of 
the  forward  sweep  the  matrices  A  +  are  averaged  between  q*j-i  and  q ^  rather  than 
homogeneously  at  the  old  iteration  level  n.  Similarly,  a  companion  backward  space 
marching  sweep  that  is  symmetric  to  equation  (5)  and  that  is  intimately  related  to 
the  backward  sweep  of  the  alternating  bidiagonal  algorithm  of  reference  (16)  is 


(/  +  1+  -  A-)6gj  =  ~A  +  Aq*)3-i  +  A~q *  -  A~  qU  +  J  (6) 


The  method  given  by  equations  (5)  and  (6)  is  von  Neumann  unconditionally  sta¬ 
ble  for  the  scalar  wave  equation.  The  analysis  shows  the  significance  of  DDADI 
approximate  factorization  in  rendering  both  the  forward  and  backward  sweeps  sep¬ 
arately  stable  regardless  of  eigenvalue  sign.  Consequently  as  the  local  Courant 
number  becomes  very  large,  :  be  robust  method  becomes  a  very  effective  (symmetric 
Gauss-Seidel)  relaxation  schei..<  ior  the  steady  equations,  a  fact  which  substantially 
contributes  to  the  very  fast  performance  that  will  be  demonstrated. 

At  a  right  computational  boundary  on  the  forward  sweep  we  solve  the  charac¬ 
teristic  boundary  point  approximation10 


(A+  +  A+)Si,n=A+1^1  -A*qnN 


,”+1 

N 


q x  and  at  a  left,  on  the  backward  sweep 


(7) 


(A  -  A  )6q ,  =  A  q™  -  A  qU  +  1  (8) 

Following  the  solution  of  equations  (7)  and  (8)  the  conservative  state  vector  is 
iteratively  corrected1''’  to  maintain  the  accuracy  of  prescribed  boundary  conditions 
while  not  disrupting  the  representation  of  the  computed  characteristic  variable:' 
running  to  the  boundary  from  the  interior.  Analysis  of  a  model  system  with  upwind 
differenced  scalar  equations  and  coupled  boundary  conditions  was  related  to  the 
linearized  bidiagonal  scheme10  by  Oligor  and  bombard18;  the  analysis  also  strongly 


supports  the  numerically  confirmed  robust  stability  of  the  present  nonlinear  method 
for  gasdynamics.  A  useful  result  of  reference  18  that  simplifies  the  procedure  of 
reference  16  is  that  on  the  forward  sweep  there  is  no  need  for  a  predictor  step  at 
the  left  boundary  J  =  1;  thus,  the  solution  sweep  begins  at  J  =  2.  Similarly,  the 
backward  sweep  begins  at  J  —  N  —  1. 

With  the  updating  at  each  step,  where  in  equation  (6)  Sq3-  =  <7^+1  —  Qji  ^  's 
clear  that  the  symmetric  pair  of  equations  (5)  and  (6)  serve  to  advance  the  solution 
two  pseudo  time  (iteration)  levels;  whereas,  the  linear  alternating  bidiagonal  sweep 
algorithm  of  reference  (16)  advances  the  solution  only  one  level.  To  maintain  con¬ 
servation  to  a  very  high  degree,  in  single  sweep  marching  in  supersonic  zones  we 
iterate  (at  least)  once  locally  at  each  space  marching  step.  The  local  iteration  serves 
to  make  the  eigenvectors  in  the  coefficient  matrices  consistent  with  the  advanced 
state  and,  thus,  provides  improved  accuracy  for  the  nonlinear  system.  It  appears 
effective  to  do  this  inner  iteration  everywhere,  i.e.  in  both  subsonic  and  supersonic 
regions,  as  the  number  of  global  iteration  steps  to  convergence  with  two  inner  iter¬ 
ations  has  been  found  reduced  by  a  factor  of  three  to  four.  Since  the  computational 
work  per  two  steps  is  about  the  same  for  the  single  level  and  two  level  schemes  and 
beyond  the  fact  that  one  saves  a  level  of  storage  in  the  space  marching  algorithm, 
the  question  arises:  Can  one  get  solutions  in  less  computational  work  through  faster 
convergence  with  the  nonlinear  space  marching  algorithm? 

One  Dimensional  Results 

First,  we  present  results  for  supersonic  flow  with  no  shock  in  Shubin’s  diverging 
nozzle.  In  purely  supersonic  zones,  the  experience  with  the  present  method  is  that 
the  solution  can  be  marched  accurately  in  one  global  iteration,  as  ought  to  be  the 
case.  Figure  4  shows  the  exact  solution  (in  solid  line)  and  the  computed  result 
from  the  first  forward  sweep.  It  is  evident  that  the  method  correctly  predicts 
the  solution  to  plotting  accuracy  in  one  forward  sweep.  With  subsequent  sweeps 
the  error  (the  difference  between  the  exact  and  the  computed  solution)  reduces  to 
machine  accuracy  in  less  than  three  global  iterations.  In  fact,  by  increasing  (from 
two)  the  number  of  inner  iterations  on  the  solution  procedure  at  each  space  marching 
step,  convergence  to  prescribed  accuracy  can  be  guaranteed  in  one  forward  sweep. 
This  is  also  true  of  contemporary  locally  linearized  unsplit  methods  in  supersonic 
flow. 
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With  the  globally  iterative  nonlinear  space  marching  formulation,  early  expe¬ 
rience  in  two  quasi  1-D  nozzle  problems  with  mixed  supersonic-subsonic  zones  is 
that  solutions  are  obtained  in  roughly  an  order  of  magnitude  fewer  iteration  steps 
than  had  been  required  with  the  previously  fast  pseudo  time  dependent  technique 
and  block  tridiagonal  solving. 

The  two  nozzle  problems  which  are  described  and  solved  by  Yee,  Beam  and 
Warming19  and  solved  with  the  CSCM  time  dependent  technique  in  references  (14) 
and  (15)  are  Shubin’s  diverging  nozzle  flow  and  Blottner’s  converging-diverging 
nozzle  flow.  Both  problems  involve  unmatched  overpressures  at  the  outflow  which 
result  in  internal  shock  terminated  supersonic  zones  and  subsonic  outflow.  For 
the  experiments  involving  flow  of  mixed  type  the  same  initial  data  given  by  Yee, 
Beam  and  Warming  -  a  linear  interpolation  between  inflow  and  outflow  values  for 
effectively  exact  solutions  of  the  problems  -  is  used  that  was  used  previously  with 
the  time  dependent  approach. 

For  flows  of  mixed  type,  in  Figures  5  and  6  respectively,  results  are  shown  for 
successive  forward  and  backward  sweeps  for  five  global  iteration  steps  with  Shubin’s 
and  Blottner’s  nozzle  flows.  In  both  cases,  the  exact  solution  as  given  by  Yee,  Beam 
and  Warming  is  shown  in  solid  line  and  the  present  computational  results  solved 
on  a  51  point  mesh,  in  boxes.  Blottner’s  nozzle  flow  is  shown  converged  after  10 
global  iteration  steps.  There  is  substantial  evidence  in  other  results  not  shown  that 
with  further  work  the  number  of  global  iterations  required  to  compute  flows  such 
as  Blottner’s  can  be  reduced  by  a  factor  of  two,  to  about  five. 

In  Figure  7,  we  show  a  subcritical,  i.e.  completely  subsonic,  flow  solution 
computed  in  only  two  global  iteration  steps  for  the  Blottner  nozzle  geometry  with 
different  inflow  conditions.  Here  the  exact  analytical  solution  derived  by  Vcnkata- 
pathy  is  shown  in  solid  line  and  our  computed  results  in  boxes. 

The  alternating  direction  sweeps  in  our  method  have  been  derived  directly  out 
of  theory  for  solving  the  implicit  set  of  difference  equations.  However,  one  can  see 
mechanistically,  numerically  speaking,  that  omitting  the  backward  sweep  from  the 
pair  and  globally  iterating  only  with  the  forward  sweep  equation  (5)  will  result  in 
permitting  the  influence  of  a  subsonic  outflow  boundary  (or  interior  disturbance)  to 
propagate  upstream  only  one  grid  point  per  global  iteration.  In  such  a  case,  which 
relates  to  other  global  iteration  methods  found  in  the  literature  and  that  also  sweep 
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only  in  the  main  flow  direction,  the  rate  of  convergence  is  greatly  inhibited  relative 
to  symmetric  sweeping  by  a  factor  of  order  roughly  the  number  of  grid  points  in  the 
subsonic  zone.  Mathematically,  this  inhibition  is  the  result  of  the  failure  to  include 
in  the  implicit  process  the  effect  of  the  eigenvectors  governing  upstream  influence 
but  to  treat  these  waves  explicitly  with  effective  CFL  unity. 

In  Figure  8  we  illustrate  the  progress  of  the  transient  solution  to  the  subcritical 
nozzle  problem  after  15  forward  sweeps,  with  the  backward  sweeps  omitted.  One 
can  clearly  see  that  the  wave  influence  of  the  outflow  boundary  has  progressed  only 
15  mesh  points  forward  of  the  outflow  boundary.  In  Figure  9  the  transient  solution 
is  shown  after  60  steps  which  is  beyond  one  characteristic  transit  time  (equivalent 
to  50  mesh  intervals)  for  the  upwind  wave  to  reach  the  inflow  boundary.  In  Figure 
10  we  show  the  history  of  the  RMS  error  in  the  primitive  variables.  The  solution  is 
found  to  converge  to  roughly  the  same  RMS  error  after  three  characteristic  times 
(150  steps)  as  the  solution  obtained  with  the  symmetric  alternating  sweep  sequence 
after  only  3  global  iterations. 

Blottner’s  supercritical  nozzle  problem  which  involves  subsonic  inflow  acceler¬ 
ating  through  a  sonic  point  to  a  supersonic  zone  terminated  by  a  shock  to  subsonic 
outflow  is  the  most  computationally  demanding  of  the  test  cases  and  indicates  the 
capability  for  the  method  to  compute  simply  and  consistently  over  the  subsonic 
forebody  and  base  regions  of  blunt  bodies  in  supersonic  flow.  Thus  the  need  for 
separate  time  dependent  codes  will  be  obviated  by  this  new  method. 

Finally,  in  Figures  11  and  12,  we  present  the  convergence  history  for  the  present 
nonlinear  scheme  and  the  linearized  time  dependent  scheme  for  completely  subsonic 
and  supersonic  nozzle  flows.  The  x-axis  shows  the  number  of  iterations  each  scheme 
requires  to  reduce  the  exact  error  to  five  orders  of  magnitude  for  various  Courant 
numbers.  It  is  evident  that  the  present  scheme  converges  extremely  fast  at  all  CFL 
numbers  compared  with  the  method  based  on  the  linearized  block  tridiagonal  solver. 

Two  Dimensional  Formulation 

For  two  dimensional  flow,  assuming  a  marching  coordinate  £,  inviscid  terms 

B+Vn+B' A„  (9a) 

and 

-B+Anq)k-i  -  B~A^q)k  (96) 
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are  added  to  the  left  and  right  hand  sides  respectively  of  both  the  forward  and  back¬ 
ward  sweep  equations  (5)  and  (6).  For  viscous  flow,  second  centrally  differenced, 
thin  layer  viscous  terms  are  also  added  in  the  77  direction  as  is  conventionally  prac¬ 
ticed,  e.g.  Steger20.  With  the  terms  for  the  tj  cross  marching  coordinate  direction, 
the  technique  now  becomes  an  implicit  method  of  lines.  Along  each  77  coordinate 
line,  one  can  solve  the  equations  coupled  with  a  block  tridiagonal  procedure.  Alter¬ 
natively,  a  further  DDADI  bidiagonal  approximate  factorization  can  be  employed  in 
the  77  direction  and  solved  either  linearly  as  in  reference  (16)  or  nonlinearly  as  here 
in  the  £  direction.  As  shown  in  the  quasi  2-D  numerical  experiments  of  reference 
(16),  DDADI  bidiagonal  approximate  factorization  is  stable  for  viscous  as  well  as 
inviscid  terms.  Finally  in  reference  (16)  there  is  a  relevant  discussion  of  the  reduced 
approximate  factorization  error  that  attends  using  DDADI  in  one  or  more  space 
directions. 

Two  Dimensional  Results 

We  present  results  for  a  45°  -  15°  axisymmetric  transonic  nozzle  flow  previously 
studied  experimentally  by  Cuffel,  Back  and  Massier21  and  computationally  by  Cline, 
Prozan,  Serra  and  Shelton  (all  referenced  in  (21))  and  ourselves16.  In  Figure  13  we 
show  results  after  10  steps  of  an  early  computation  run  at  a  local  CFL  number  of 
14  with  the  present  first  order  single  level  scheme.  Except  for  the  addition  of  an 
error  correction  procedure15  to  counter  numerical  inflow  boundary  condition  drift, 
a  factor  which  has  improved  the  present  solution  in  the  vicinity  of  the  axis,  the 
effectively  converged  results  found  here  are  the  same  as  those  given  for  the  two 
level  scheme  in  reference  16.  (As  long  as  the  problem  has  a  unique  solution,  the  two 
schemes  must  give  equivalent  results  since  the  right  hand  side  difference  equation 
sets,  including  boundary  approximations,  are  the  same.) 

For  the  solution  given  in  Figure  13,  we  noted  a  very  rapid  rate  of  reduction 
in  residual,  three  orders  of  magnitude  in  ten  steps.  This  compares  with  60  steps 
given  in  reference  (16)  for  the  solution  obtained  with  the  two  level  scheme.  The 
rapid  convergence  found  in  this  transonic  problem  for  the  CSCM-S  method  with 
viscous  terms  provides  the  reasonable  expectation  of  similar  fast  results  to  be  ob¬ 
tained  without  viscous  effects.  Thus  the  method  in  multidimcnsions  appears  to  have 
attractive  potential  for  an  improved  transonic  Euler  solver  as  well  as  Navier-Stokes 
solver. 


Next,  we  present  first  order  inviscid  and  viscous  results  for  an  inlet  problem 
shown  in  Figure  14.  The  pressure  contours  for  the  first  order  inviscid  solutions 
are  shown  in  Figure  15.  Figure  16  shows  the  first  order  viscous  results.  The  vis¬ 
cous  computation  shows  the  presence  of  the  leading  edge  shock.  The  flow  structure 
compares  very  well  with  the  theoretical  (for  the  inviscid  case)  and  other  computa¬ 
tional  results.  In  Figure  17,  the  inviscid  and  viscous  wall  pressure  are  compared 
with  the  exact  solution  (inviscid).  Figure  18  shows  the  convergence  history  of  the 
RMS  residue  of  all  the  conservative  flow  variables  for  the  inviscid  problem  solved 
at  CFL  =  100  with  4  inner  iterations  at  every  axial  location.  For  the  inviscid  case, 
only  forward  marching  was  carried  out  and  backward  marching  was  omitted.  The 
solution  has  converged  for  practical  purposes  at  the  end  of  the  first  sweep.  The 
residue  reaches  machine  accuracy  in  10  iterations.  In  reference  22,  we  show  the 
residual  reduction  versus  inner  iteration  number  in  single  sweep  solutions  for  super¬ 
sonic  flow  and  compare  results  with  contemporary  PNS  procedures.  Finally  in  work 
described  in  reference  23  we  have  applied  the  single  level  scheme  to  the  solution  of 
the  coupled  forebody  captured  shock  layer  and  massively  separated  wake  flow  of  a 
hypersonic  axisymmetric  AOTV. 

Three  Dimensional  Method  of  Planes  Algorithm 

In  reference  24  we  presented  a  symbolic  algebra  for  DDADI  approximate  fac¬ 
torization  and  derived  single  level  relaxation  schemes.  The  algebra  is  based  on  the 
implicit  difference  stencil  of  the  implicit  method.  Here  we  will  show  how  the  ap¬ 
proach  can  be  used  to  derive  an  symmetric  Gauss-Seidel  implicit  method  of  planes 
relaxation  algorithm. 

The  unfactored  three  dimensional  linearized  implicit  method  can  be  represented 
by  the  symbolic  matrix  expression 

B~ 

I  ^  C~ 

-A+  -  D - A~ 

-C+  "  I 
- B+ 

The  block  diagonal  where  matrix  D  —  I  +  A+  —  A~  +  B+  —  B~  +  C+  —  C~ .  A 
once  DDADI  approximate  factorization  in  the  f  coordinate  direction  leads  to  the 
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expression 


-C+' 


-B+ 


-B+ 


By  analogy  with  the  derivation  of  the  single  level  scheme  for  quasi  1-D  flow  from 
DDADI  bidiagonal  approximate  factorization  we  identify  the  above  expression  with 
the  alternate  space  marching  implicit  method  of  planes  algorithm 

Forward  Sweep 


-c+  d  C-  =  rhs\  ?”+,■  , 

-s+ 

Backward  Sweep 


-C+  D  C~  6q^+1  =  RHS[  q»+1  ,  <£+2  ] 

~B+  J 

In  the  planes  the  coupled  block  matrix  problem  can  be  further  simplified  by  the 
approximate  factorization 


[  -C+  ,  D  ,  C~  ]  D~l  [  ~B+  D  B~  ]  =  RHS 

which  leads  directly  to  the  block  tridiagonal  solution  sequence 

(  -C+  D  C~  ]<V  =  RHS 
[  -B~  D  B~  \6g  =  D8q * 

Three  Dimensional  Results 

This  3-D  space  marching  algorithm  has  been  tested  against  an  axisymmetric 
viscous  flow  problem  of  a  supersonic  RL-10  nozzle  exhausting  over  a  backward  facing 
step  into  a  cylindrical  shroud. 

In  Figure  10  we  show  a  3-D  perspective  view  of  the  wall  surface  mesh  in  the 
quarter  sector  for  which  we  solved.  A  section  of  computational  mesh  in  a  longitu¬ 
dinal  plane  through  the  axis  is  shown  in  Figure  20.  For  such  planes,  Figures  21  and 


22  show  pressure  contour  and  velocity  vector  plots  that  exhibit  the  expected  flow 
structures  -  a  weak  shock  in  the  nozzle  diffuser,  a  strong  expansion  at  the  nozzle 
exit  with  subsequent  recompression  shock  off  the  shroud  wall  and,  lastly,  a  substan¬ 
tial  region  of  streamwise  separated  flow  under  the  backward  facing  step.  The  Mach 
contour  plot  of  Figure  23  for  a  cross  flow  data  plane  in  the  shroud  shows  the  ex¬ 
cellence  with  which  the  method  and,  in  particular  the  lateral  boundary  conditions, 
numerically  maintain  the  axisymmetry  of  the  solution. 

Patched  Mesh  Boundary  Procedures 

In  previous  work  Lombard,  et  al15,16  and  Oliger  and  Lombard18  gave  stable 
implicit  procedures  for  computing  the  solution  at  external  boundaries  of  a  com¬ 
putational  domain.  Those  procedures  generalized  the  work  of  Kenzer  to  matrix 
coupled  linearized  boundary  conditions  complementing  a  set  of  advective  difference 
equations  (associated  with  well  posed  characteristics)  to  the  boundary. 

Under  the  present  contract  we  have  explored  the  problem  of  implicitly  cou¬ 
pling  at  interior  patch  boundaries  the  global  solution  on  a  system  of  multiple  patch 
meshes.  The  approach  we  have  taken  in  this  research  is  numerical  experimenta¬ 
tion  among  a  number  of  boundary  treatment  approximations  to  the  solution  of  the 
continuous  domain  problem.  For  comparison  purposes  with  previous  single  mesh 
results,  numerical  experiments  were  performed  with  the  well  tested  two  level  pseudo 
time  relaxation  CSCM  scheme  at  the  interior  points  of  the  mesh  patches. 

Quasi  2-D  Studies 

Our  first  experiments,  to  be  described  here,  dealt  with  breaking  a  single  co¬ 
ordinate  line  into  segments  and  solving  sequentially  on  each  the  equations  for  a 
quasi  2-D  viscous  compressible  flow.  The  model  problem,  with  which  we  experi¬ 
mented  in  reference  16  for  an  uninterrupted  mesh,  is  a  transient  pipe  flow  resulting 
from  an  initially  nonequilibrated  pressure  between  the  axis  and  wall  boundaries. 
Three  kinds  of  cases  were  run  with  the  two-level  linearized  implicit  procedure;  all 
featured  sequential  solving  on  patches  with  at  least  one  point  of  overlap.  Case  1 
had  frozen  boundary  data,  obtained  from  the  solution  on  neighboring  meshes  in  an 
operationally  explicit  manner.  Specifically  at  left  and  right  first  computed  interior 
points  of  a  patch,  we  solved  the  bidiagonal  equations  respectively 

(/  +  -  A~A)6q2  =  A+qni  -  (A+  -  A~)q”  -  A'q*  (10a) 
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and 


[I -A  +  A  +  V)<^_i  =  A+qNU_  2  -  {A+  -  A  )qNn_  1  ~  A  (106) 

Here  the  symbols  n  indicate  the  “frozen”  boundary  data  coming  from  the  solution 
on  the  interior  of  an  adjoining  and  partially  overlying  patch  was  at  the  old  iteration 
level  in  the  global  procedure.  In  Case  1  the  solution  on  all  the  grids  was  effectively 
updated  at  the  same  time. 

Case  2  featured  reverse  sequential  cycling  of  the  two  level  time  dependent  solu¬ 
tion  procedure  through  the  grids  on  alternate  global  steps.  In  the  forward  sequence 
the  right  interior  patch  boundaries  were  (a)  frozen  or  (b)  computed  with  one-sided 
characteristic  boundary  conditions  obviating  any  change  in  the  characteristic  data 
from  outside  (to  the  right  of)  a  patch.  The  left  interior  patch  boundaries  inherited 
implicit  data  ( 6q )  from  the  solution  on  the  computed  patch  to  the  left.  In  the 
backward  sequence  all  the  roles  were  reversed  including  the  directions  of  forward 
elimination  and  back  substitution  in  the  tridiagonal  matrix  inversion  procedure. 

Case  3  featured  cycling  through  the  patches  in  the  predictor  (forward  elimi¬ 
nation)  step  of  the  solution  procedure,  inheriting  implicit  left  boundary  data  as  in 
Case  2.  Then  the  patches  were  cycled  through  in  reverse  order  on  the  corrector 
(back  substitution)  step.  The  result  (save  for  interpolation  if  data  points  of  the 
grids  were  interlaced)  is  identically  equivalent  to  solving  on  an  uninterrupted  single 
patch  mesh. 

The  results  of  the  three  sets  of  experiments  were  comparable  for  local  time 
steps  based  on  constant  CFL  number  up  to  about  50.  Beyond  that  the  rate  of 
convergence  of  the  solutions  for  Cases  1  and  2  diminished  and  at  sufficiently  high 
CFL  number  failed  to  converge.  The  efTort  involved  in  Case  2  with  the  two  level 
scheme  was  not  rewarded  relative  to  the  simple  procedure  of  Case  1.  Within  the 
framework  of  the  single-level  symmetric  Gauss-Seidel  implicit  relaxation  scheme, 
which  is  linearly  equivalent  to  the  bidiagonal  scheme  employed  at  the  boundary, 
Case  2a  is  operationally  no  more  difficult  than  Case  1  and  more  closely  approximates 
Case  3.  Case  2b  has  a  consistency  problem  that  inhibits  firm  convergence. 

Figure  24  compares  rms  residual  (density)  convergence  history  for  Cases  1  and 
3  with  three  mesh  segments  with  two  cells  of  overlap  and  run  at  CFL  25.  As  one 
might  expect,  the  effectively  uninterrupted  mesh  procedure  is  found  to  converge 
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(two  to  three  times)  faster  for  about  the  first  50  steps  through  the  major  transient. 
After  that  the  performance  is  comparable.  The  performance  of  the  frozen  boundary 
treatment  Case  1  with  five  segments  at  CFL  25  was  found  to  be  not  materially  worse 
than  with  three  segments,  but  the  performance  degradation  with  increasing  CFL 
was  found  to  proceed  faster  with  increasing  number  of  patches. 

From  the  results  of  the  quasi  2-D  experiments  and  extrapolation  from  experi¬ 
ence  with  the  single  level  scheme,  we  observe  that  the  simple  boundary  procedure 
with  frozen  conservative  variable  data  taken  from  the  solutions  at  either  iteration 
level  n  or  n  -f  1  on  adjacent  meshes  and  coupled  implicitly  by  alternating  direction 
sequential  solving  through  the  patches  has  robust  stability  to  sufficiently  high  CFL 
number  to  yield  a  rate  of  convergence  meeting  our  needs. 

Two  Dimensional  Patch  Boundary  Treatment 

In  earlier  work  solving  the  Euler  equations  on  multiple  patch  meshes  Benek, 
et  al.25  employed  linear  interpolation  of  the  conservative  variables  from  the  interior 
solution  of  one  mesh  to  give  Dirichlet  boundary  data  for  the  other.  With  several 
points  of  mesh  overlap  at  the  mesh  boundaries,  transonic  solutions  obtained  with 
central  differencing  exhibited  considerable  oscillation  in  the  vicinity  of  shocks  prop¬ 
agating  through  the  boundary  region.  Eberhardt26  with  the  code  of  Benek,  et  al.25 
attempted  to  reduce  the  oscillation  and  attendant  stability  problems  encountered 
in  the  vicinity  of  a  bluff  body  shock  intersecting  an  embedded  patch  boundary  by 
introducing  a  characteristic  computed  boundary  point  approximation  with  scalar 
upwind  difference  equations  in  Riemann  variables.  In  his  procedure,  interpolation 
was  performed  in  only  the  variables  whose  characteristics  ran  to  the  patch  boundary 
from  outside  the  computational  domain.  Eberhardt  based  the  decision  about  do¬ 
main  of  dependence  of  the  characteristics  on  eigenvalues  computed  within  the  patch 
domain.  When  this  decision  was  compatible  with  the  flow,  then  the  characteris¬ 
tics  boundary  procedure  gave  markedly  superior  results  compared  to  interpolation 
of  conservative  variable  data,  which  leads  to  solution  overspecification  in  subsonic 
zones.  In  other  cases  where  incorrect  domain  of  dependence  was  inferred,  the  char¬ 
acteristic  boundary  procedure  was  unstable. 

Here  based  on  our  quasi  2-D  studies  described  above,  we  give  a  simple,  robustly 
stable  implicit  approach  to  computing  solutions  of  the  conservative  equations  of 
gasdynamics  on  either  composite  or  overset  meshes.  Without  requiring  special  flux 
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conservative  operators,  but  rather,  interpolating  conservative  variable  data  at  mesh 
boundaries,  the  implicit  upwind  method  is  accurate  and  relatively  free  of  oscillation 
where  shocks  intersect  interior  patch  boundaries.  In  a  supersonic  inlet  problem 
with  expansion  and  reflected  shocks,  we  exhibit  the  capability  to  conveniently  carry 
out  with  rapidly  convergent  implicit  methods  for  systems  of  equations  the  adaptive 
refined  meshing  strategy  in  overset  patches  proposed  by  Berger  and  Oliger2.  Fur¬ 
ther,  the  test  case  shows  concretely  in  a  realistic  aerodynamic  problem  the  savings 
in  mesh  points  (about  an  order  of  magnitude  here  in  two  dimensions)  for  similar 
accuracy  that  flow  structure  aligned  adaptive  patched  meshing  affords  compared  to 
uniform  grid  refinement. 

The  factors  in  our  approach  are  supported  by  previous  research  by  ourselves  and 
generically  by  others  cited  in  reference  27  and  are  proven  in  numerical  experiments 
reported  here  and  elsewhere28.  We  employ  an  implicit  conservative  upwind  scheme 
CSCM1G  with  which  in  the  present  work  we  can  solve  to  either  first  or  second 
order  spatial  accuracy  the  Euler  or  compressible  Navier-Stokes  equations  in  two- 
dimensional  planar  or  axisymmetric  flows.  The  flux  difference  split  upwind  schemes 
of  generalized  Roe  form  such  as  CSCM  have  a  number  of  qualities  that  make  them 
ideal  for  the  purpose  of  solving  on  discontinuous  patched  mesh  systems. 

First,  conservative  schemes  in  the  Roe  form  satisfy  Roe’s  property  U  that 
guarantees  the  correct  speed  for  captured  shocks.  The  CSCM  scheme  has  been 
tested  in  a  wide  variety  of  internal  and  external  transonic  and  hypersonic  flows16,22 
and  has  been  found  to  capture  strong  and  weak  shocks  accurately  in  location  and 
with  little  oscillation.  The  shock  transition  is  particularly  sharp,  about  two  mesh 
cells  wide,  on  an  aligned  grid;  and  this  factor  will  be  accommodated  as  much  as 
possible  in  our  adaptive  patched  mesh  strategy. 

Second,  in  the  Roe  form  the  difference  operators  on  conservative  variable  data 
represent  the  effects  of  differencing  to  characteristic  data  only  for  disturbances  prop¬ 
agating  toward  the  given  node  and  reject  the  mathematically  unstable  contribution 
from  disturbances  that  may  be  propagating  downwind  of  the  node.  The  one  sided 
upwind  difference  operation  represents  identically29  the  (split)  conservative  partial 
flux  difference  between  the  nodes.  To  the  extent  that  the  data  from  an  adjacent 
mesh  is  consistent  with  the  solution,  then  the  associated  upwind  partial  flux  differ¬ 
ence  to  that  data  will  serve  to  provide  at  convergence  consistency  of  the  partial  flux 
convective  into  the  given  mesh  from  its  neighbor,  and  vice  versa  for  signals  of  the 


other  eigenvalue  sign.  Thus  the  method  acts  within  truncation  error  to  provide  the 
similar  continuity  of  the  flux  tensor  among  patched  grids  that  exists  in  the  phys¬ 
ical  solution  across  shocks  and  would  obtain  on  a  single  grid  alone.  The  correct 
domain  of  dependence  coupled  with  the  well  posed  characteristic  boundary  point 
approximations16  tend  not  to  support  oscillatory  disturbances  but  convect  them 
out  of  the  flow  domain. 

Third,  one  sided  difference  interval  averaged  eigenvalues  let  majority  rule  de¬ 
termine  the  direction  of  local  signal  propagation.  When  applied,  as  we  do,  to  a 
difference  operator  between  boundary  data  obtained  from  an  adjacent  mesh  and 
the  local  mesh  solution  point,  the  data  of  both  meshes  participate  in  making  the 
decision  as  to  whether  an  incoming  signal  is  being  sent.  Both  in  concept  and  our 
experience,  this  factor  seems  to  overcome  the  inter  mesh  communication  difficulty 
experienced  by  Eberhardt26  with  his  characteristic  boundary  procedure. 

Lastly,  with  the  CSCM  difference  equations,  with  diagonally  dominant  approx¬ 
imate  factorization16,22  that  retains  on  the  diagonal  the  contributions  from  both 
sets  of  eigenvalues  in  what  is  effectively  an  absolute  value  Jacobian  matrix,  we  can 
solve  the  equations  either  with  two  data  level  linearized  block  implicit  methods16 
or  with  a  single  data  level  relaxation  technique22  that  is  substantially  more  rapidly 
convergent  than  the  linearized  implicit  procedures.  As  can  be  inferred  from  the 
theory  and  numerical  experiments  of  reference  22,  the  use  of  DDADI  on  the  solu¬ 
tion  point  while  differencing  effectively  explicitly  to  data  obtained  from  an  adjoining 
mesh  (which  may  be  at  either  the  old  (n)  or  new  (n+1)  level)  is  unconditionally  sta¬ 
ble.  When  the  solutions  on  the  patched  meshes  are  alternately  updated  using  either 
the  linearized  implicit  or  relaxation  methods  the  global  procedure  is  implicit.  As 
will  be  shown  here,  the  robust  stability  of  the  global  procedure  has  been  confirmed 
to  approximately  coarse  mesh  CFL  100. 

One  point  that  has  not  been  touched  on  is  the  form  of  interpolation  that  we 
use.  Differencing  to  interpolated  data  is  equivalent  to  a  weighted  sum  of  differences 
operating  on  the  interpolants.  It  is  intuitive  that  for  robust  stability  each  of  these 
assumed  upwind  differences  ought  to  be  well  posed.  This  implies  that  the  inter¬ 
polation  weights  should  all  be  positive  and  the  domains  of  dependence  of  all  the 
interpolants  should  be  outside  (i.e.  on  the  assumed  side)  of  the  solution  point  with 
respect  to  its  mesh  interior.  Neither  of  these  properties  was  shared  by  the  data 
interpolation  schemes  used  by  Benek,  et  al.25  or  Eberhardt26. 
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Finally,  extending  a  direction  of  Benek,  et  al.25,  we  do  not  compute  on  sections 
of  coarser  meshes  underlying  patches  of  overset  refinement.  Our  data  structure  and 
automated  procedures  for  the  consequent  partitioning  of  meshes  are  described  in 
reference  30.  The  partitioning  concept  in  which  coordinate  lines  of  a  patch  have 
differing  (index)  lengths  in  computational  space  also  leads  to  useful  possibilities 
including  as  will  be  shown  here,  fitting  mesh  patches  obliquely  to  boundaries,  e.g. 
to  sharply  capture  reflecting  shocks. 

Interior  Boundary  Treatment  -  First  Computed  Point  Formulation 

We  review  here  the  first  order  scheme  for  one  space  dimension.  In  equation  (2) 
the  CSCM  flux  difference  splitting  is 

( A+  +  A~)Aq  =  A  F+  +  A  F~  =  A  F 

with 

A±  =  (MT/±T_1M_1)I=  A±A  (11) 

and 

=  ^(/  ±  sgn(A)) 

exhibiting  the  similarity  transformation  that  diagonalizes  the  constructed  flux  dif¬ 
ference  Jacobian  A.  Here  A  is  a  diagonal  matrix  of  the  interval  averaged  eigenvalues 
that  through  the  truth  function  diagonal  matrices  /±  make  the  decisions  about  di¬ 
rections  of  characteristic  wave  propagation  and  whether  or  not  to  send  signal  to 
the  solution  point.  Thus  in  the  equation,  A  +  Ag)y_i  represents  the  convection 
of  characteristic  wave  contributions  in  the  positive  coordinate  direction  from  grid 
point  j  —  1  to  solution  point  j  and  A~ ,  in  the  negative  direction  from  j  +  1  to 
j.  As  the  result  of  incorporating  multiplicatively  the  (local)  time  step  (for  pseudo 
time  relaxation)  and  the  spatial  (divided)  differences  in  the  matrices,  the  numerical 
eigenvalues  are  Courant  numbers  for  the  characteristic  waves  whose  speeds  are  u, 
u  +  c  and  u  -  c,  with  c  the  sound  speed. 

For  exterior  boundaries  of  the  computational  domain,  at  a  right  boundary  point 
we  retain  only  the  right  running  positive  eigenvalue  contributions  and,  negative 
at  left.  At  such  boundary  points  the  identity  matrix  on  the  left  hand  side  of  the 
equation  (1)  is  replaced  by  a  matrix10  that  contains  the  time  linearized  conservative 
variable  representations  of  the  computed  characteristic  variable  fluctuations  at  the 
boundary  and  also  imposed  boundary  condition  relations. 


Of  more  interest  here  is  what  we  do  at  an  interior  patch  boundary.  The  left 
hand  side  of  equation  (2)  has  the  tridiagonal  structure 


I  +  A+  -  A- 


In  relation  (12),  the  central  block  which  we  call  D  can  be  seen  to  contain  the 
absolute  values  of  the  eigenvalues  for  signals  propagating  to  the  solution  point  from 
either  left  or  right.  Indeed,  the  simplified  approximation  to  equation  (2) 


nr  r+  n,n  +  1 
Doq-i  —  A  q 

3  3-  1 


r_,  n  ~  n,n  +  1 
A  <7  .  -  A  q 

3  3  +  1 


leads  to  an  operationally  explicit  implicit  relaxation  procedure22  that  is  uncondi¬ 
tionally  stable  either  as  a  computed  interior  patch  boundary  point  or  general  interior 
point  relation.  Here  n,  n  +  1  means  data  from  either  time  level.  If  the  interior  point 
implicit  solution  procedure  is  two  level,  then,  the  term  of  equation  (12)  at  the  inte¬ 
rior  point  j  —  1  or  j  -\- 1  will  be  linearized  (assumed  at  the  n  +  1  level)  as  in  equations 
(10a)  or  (10b). 

For  left  or  right  boundary  points,  the  frozen  (i.e.  not  computed  on  the  patch) 
data  at  j  —  1  or  j  +  1  in  equation  (13)  is  gotten  from  adjacent  patch  mesh  data. 
If  the  mesh  system  is  composite  and  mesh  lines  cross  the  boundary  to  the  solution 
point,  then,  the  frozen  boundary  data  is  the  solution  point  data  of  the  adjoining 
mesh.  For  the  case  of  nodes  on  lines  ending  at  the  patch  boundary,  which  case 
relates  equally  to  composite  mesh  with  lines  of  either  patch  ending  at  the  boundary 
or  to  overset  meshes,  the  frozen  data  is  got  by  interpolation  of  adjacent  mesh  patch 
data  to  the  boundary  point  location. 

Linear  interpolation  to  an  included  point  on  a  coordinate  line  or  within  a 
polygonal  cell  involve  only  positive  weights  on  the  interpolant  data.  In  most  of 
the  numerical  experiments  made  to  date  with  overset  grids  we  have  employed  a 
bilinear  interpolation30  based  on  the  four  corner  points  of  the  overlapping  mesh  cell 
enclosing  the  frozen  boundary  point.  However,  with  less  data  processing  a  linear 
interpolation  involving  the  three  corner  points  on  the  including  triangle  (Figure 
25a)  generalizes  to  the  use  of  the  four  corner  points  of  the  enclosing  tetrahedron  in 
three  dimensions. 


In  the  composite  mesh  case,  Figure  25b,  interpolation  is  naturally  along  an 
interior  coordinate  line  paralleling  the  zonal  boundary.  Such  interpolation  is 
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one-dimensional  for  a  two-dimensional  problem  and  two-dimensional  for  a  three- 
dimensional  problem.  The  generic  composite  grid  problem  just  described  and  which 
we  have  tested  among  the  numerical  experiments  to  be  reported  in  the  next  section 
also  serves  as  a  gedanken  model  problem  for  the  overset  mesh  case.  Possible  solu¬ 
tions  that  come  to  mind  by  analogy  are  shown  in  Figure  25c  and  Figure  25d.  In 
both  figures  the  interpolation  is  along  two  point  coordinate  lines  segments  (or  three 
point  triangular  surfaces  in  3-D)  of  the  adjacent  mesh  and  thus  is  a  direct  analog 
with  the  attendant  data  requirements  of  the  composite  mesh  case. 

We  use  stable  and  consistent  first  order  differencing  and  interpolation  at  the 
boundaries  regardless  of  the  order  of  accuracy  of  the  difference  approximation  in  the 
patch  interiors.  Since  the  divided  differences  of  the  computed  boundary  point  ap¬ 
proximation  are  of  the  same  accuracy  as  the  interpolation,  the  approaches  sketched 
in  Figures  25c  and  25d  may  be  regarded  as  letting  the  difference  operator  perform 
the  interpolation  (to  uniformly  spaced  data  Figure  25a)  in  the  direction  away  from 
the  computed  boundary.  Thus,  to  the  extent  that  the  solution  is  locally  well  rep¬ 
resented  by  a  linear  function  the  approaches  sketched  in  Figures  25a  and  25c  and 
25d  are  equivalent.  The  treatment  shown  in  Figure  25a,  however,  requires  the  same 
dimensionality  of  the  (triangular)  interpolation  procedure  as  that  of  the  problem 
and  one  higher  than  for  the  (linear)  treatment  of  Figures  25c  and  25d. 

As  a  final  theoretical  point  regarding  data  exchange  at  patch  mesh  boundaries, 
we  note  here  that  equations  (ll)  imply20 

A  F±  =  A±AF  (14) 

Hence  we  may  equivalently  write  the  right  hand  side  difference  operators  of  equation 
(13)  directly  in  terms  of  flux  components.  An  advantage  of  this  approach  is  that 
the  flux  components  normal  to  the  shock  discontinuity  are  continuous  across  the 
discontinuity.  Thus  interpolating  flux  data  from  one  grid  to  serve  as  needed  bound¬ 
ary  data  of  another  can  be  a  smoother  more  accurate  procedure  than  interpolating 
conservative  variable  data. 

We  close  this  section  by  noting  that,  consistent  with  the  2-D  interior  point 
schemes,10'22  we  difference  along  the  computed  boundary  coordinate  line  with  op¬ 
erators  written  for  the  boundary  aligned  coordinate.  In  the  diagonally  dominant 
approximate  factorizations  that  we  employ  in  multidimcnsions,  the  convective  ma- 


trices  with  absolute  values  of  the  eigenvalues  for  both  coordinate  directions  are 
retained  in  the  diagonal  block. 

Two  Dimensional  Numerical  Experiments 

To  achieve  both  accuracy  and  robust  stability  in  this  difficult  problem  area, 
particularly  with  overset  meshes,  theory  can  provide  insight  into  what  to  attempt 
but  the  acid  test  of  numerical  methods  is  performance  in  relevant  numerical  exper¬ 
iments. 

We  present  here  some  sample  results  with  major  findings  from  a  substantial 
number  of  experiments27  designed  to  test  various  aspects  of  the  accuracy  and  sta¬ 
bility  question  for  the  conservative  system  of  equations  for  gasdynamics  solved  on 
patched  meshes. 

The  numerical  experiments  to  be  discussed  here  involve  solution  of  an  inviscid 
flow  in  a  Mach  5  inlet  with  10°  compression  ramp  that  we  have  employed  in  previous 
experiments  with  first  and  second  order  upwind  methods  on  uniform  meshes10.  The 
problem  involves  two  of  the  generic  kinds  of  flow  structure,  shock  and  expansion 
fan,  which  are  not  possible  to  resolve  both  efficiently  and  to  the  extent  desired  on 
uniform  mesh.  As  the  result  of  interaction  of  the  expansion  fan  with  the  compression 
corner  and  reflected  shocks,  they  curve  in  non  simple  regions  for  which  the  exact 
solution  is  not  known  analytically.  The  coarse  base  level  grid  of  the  experiments 
has  26  x  26  points. 

Composite  Grid  with  Boundary  Overlap 

As  a  simple  test  of  employing  frozen  interpolated  boundary  data,  we  show  in 
Figure  26a  and  26b  the  grid  and  density  contours  for  a  patched  mesh  with  two 
full  cell  overlap  and  refinement  with  twice  as  many  mesh  points  in  the  streamwise 
direction  in  the  lower  patch.  Thus  every  other  mesh  point  at  the  upper  interior 
boundary  of  the  lower  patch  is  interpolated  between  computed  conservative  variable 
data  of  the  upper  patch  along  the  common  streamwise  coordinate  line.  This  test 
case  solving  sequentially  on  the  patches  with  the  first  order  method  is  numerically 
stable  with  local  time  steps  based  on  CFL  100.  There  is  no  oscillation  in  the  solution 
in  the  vicinity  of  the  patch  boundary  (shown)  and  refinement  in  the  lower  patch 
has  served  to  sharpen  the  solution  in  that  region,  though  high  gradient  regions  of 
the  solution  are  very  smeared  on  the  coarse,  nonaligncd  meshes. 


Solutions  on  Uniformly  Refined  Mesh 

As  a  standard  for  comparison  with  results  from  overset  patched  refinement,  we 
show  in  Figures  27a  and  27b  pressure  contours  for  first  and  second  order  upwind 
methods  on  101  x  101  point  uniform  grids,  i.e.  4  times  refined  in  each  coordinate 
direction. 

Adaptive  Refinement  in  Overset  Grids 

Based  on  a  uniform  coarse  mesh  solution  similar  to  Figure  26b,  the  coarse  mesh 
was  overlaid  with  two  refined  mesh  patches,  Figure  28a,  aligned  with  the  compres¬ 
sion  corner  and  reflected  shocks.  Note  in  the  reflection  r>  O.on,  the  two  overset 
patches  have  been  constructed  to  share  the  same  coordinate  lines  for  superior  grid 
communication.  The  coarse  grid  is  segmented  (broken)  under  the  overset  patches 
and  the  refinement  is  segmented  to  terminate  at  the  reflecting  boundary  (symmetry 
plane).  In  Figures  28b  and  28c  we  show  respectively  pressure  contours  for  the  first 
and  second  order  upwind  methods  on  the  overset  grid. 

Discussion 

While  the  adaptive  refinement  about  doubled  the  coarse  mesh,  the  shock  struc¬ 
tures  treated  are  better  resolved  than  with  the  uniformly  refined  mesh  in  16  times 
the  points.  Thus  the  results  demonstrate  an  order  of  magnitude  improvement  in 
data  efficiency  to  be  gained  by  overset  refinement.  We  sketch  in  Figure  28d  a 
strategy  of  refinement  for  the  as  yet  untreated  expansions  of  the  problem. 

Graphics  for  Patched  Grids 

Graphics  is  an  important  tool  to  develop  and  debug  numerical  codes  and  an¬ 
alyze  numerical  results.  In  some  problems,  contour  plots  of  certain  flow  quantities 
such  as  pressure,  Mach  number  are  sufficient  to  look  at;  and,  in  other  problems,  one 
may  need  to  look  at  velocity  vectors  to  understand  the  solution  better.  In  the  case 
of  multiple  grids,  complication  arises  due  to  the  fact  that  solutions  in  more  than 
one  grid  are  available  in  some  regions  of  the  flow  domain  and  the  global  solution 
needs  to  be  composed  from  solutions  on  individual  grids. 

Graphics  for  Single  Grid  Solution 

Before  we  deal  with  multiple  grid  graphics,  we  outline  how  we  analyze  the 
single  grid  solutions.  Since  multiple  grid  solutions  are  made  up  of  single  grid  solu¬ 
tions,  insight  into  the  single  grid  graphics  will  help  to  understand  the  multiple  grid 
graphics  plotting  strategy.  In  the  process  of  solving  problems  there  are  three  roles 
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for  graphics.  First,  the  grid  and  the  starting  solution  can  be  checked.  Second,  as 
the  solution  evolves,  graphics  can  be  used  as  a  tool  to  debug  and  understand  the 
time  evolution  of  the  solution  process.  Third,  the  converged  or  the  final  solution  is 
plotted  to  study  the  physics  and  compare  the  present  solution  with  other  solutions. 
Most  of  the  above  tasks  can  be  grouped  into  three  categories  of  plotting:  1)  Grid 
plotting;  2)  Contour  plotting  of  various  scalar  flow  quantities  and  velocity  vector 
plotting;  3)  Comparison  plotting.  Here  we  deal  with  the  first  two. 

Grid  Plotting 

Single  grid  plotting  is  done  by  essentially  drawing  straight  line  segments.  We 
use  the  commercial  macro  plotting  package  DISSPLA*  for  all  plotting  purposes. 
Our  plot  program  links  the  macro  DISSPLA  program.  Only  higher  level  commands 
need  to  be  defined  on  the  plotting  program  and  all  the  lower  level  commands  are 
defined  in  the  DISSPLA  program  and  are  transparent  to  the  users.  Given  the 
ordered  set  of  grid  points,  the  grid  plotting  routine  calls  the  line  segment  drawing 
command  for  every  line  segment  and  plots  the  grid.  Since  only  discrete  grid  values 
are  used  in  all  the  finite  difference  calculations,  the  representation  of  the  grid  by 
linear  elements  is  most  appropriate  and  no  interpolation  or  smoothing  is  necessary. 

The  same  philosophy  is  adhered  to  in  plotting  contour  curves  of  flow  quantities. 
Independent  of  the  order  of  the  scheme  and  the  solution  accuracy,  the  solution  is 
available  only  at  discrete  node  points  and  any  attempt  to  represent  the  solution  in 
a  smoother  sense  can  only  corrupt  the  solution.  With  this  in  mind,  we  use  a  fast, 
simple  and  robust  contour  plotting  routine. 

Contour  Plotting 

Contour  plots  for  the  two  dimensional  problems  are  very  useful  to  show  how 
accurately  shock  and  other  flow  structures  have  been  captured.  Any  scalar  flow 
quantity,  such  as  pressure,  density,  Mach  number,  temperature,  stream  function  or 
vorticity  can  be  plotted  for  various  constant  contour  levels.  The  general  method 
for  plotting  contour  levels  are  given  in  the  following  section.  First,  the  desired  flow 
quantities  are  calculated  from  the  dependent  variable  solution  vector  at  all  grid 
points  and  the  contour  subroutine  is  called  with  the  set  of  contour  values.  The 
contour  subroutine  computes  and  plots  the  various  contour  curves,  usually,  in  the 
physical  domain. 

‘DISSPLA  is  a  trademark  of  Integrated  Software  Systems  Corp. 
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Since  any  finite  difference/finite  volume  formulation  solves  the  flow  field  in  an 
ordered  set  of  grid  points/cells,  the  construction  and  the  execution  of  the  contour 
program  is  structured  on  the  basic  grid  cell.  We  draw  all  contour  curves  cell  by  cell 
as  we  sweep  through  the  complete  grid.  At  the  corners  of  a  cell,  flow  function  F 
has  values  FI,  F 2,  F 3,  F4.  It  is  desired  that  the  contour  curve  for  function  level 
FC  need  to  be  plotted.  Then  along  each  side  of  the  base  grid  cell,  where  cross  over 
points  of  the  contour  curve  FC  occur  they  can  be  found  through  linear  interpolation. 
By  connecting  sequentially  the  cross  over  points  encountered  among  the  sides,  one 
obtains  the  part  of  the  contour  curve  or  curves  in  a  given  cell  correpsonding  to 
the  contour  value  FC.  By  repeating  this  process  for  all  cells,  the  complete  set  of 
contour  curves  for  the  whole  domain  can  be  found. 

Since  the  interpolation  used  to  find  the  contour  curve  in  the  base  cell  is  only 
linear,  the  accuracy  of  the  contour  curve  inside  the  cells  are  also  linear.  This  does 
not  mean  that  the  solution  represented  by  the  contour  curves  is  first  order.  The  dis¬ 
crete  solution  at  the  nodes  is  as  high  an  order  as  the  solution  technique.  If  one  used 
higher  order  interpolation  to  represent  the  discrete  data,  then  the  contour  curves 
represent  not  just  the  solution  but  non-physical/extra  smoothing.  Also  higher  or¬ 
der  interpolations  to  represent  discrete  data  can  result  in  an  oscillatory  solution 
representation  when  such  is  not  the  case. 

Velocity  Vector  Plotting 

Apart  from  the  contours,  at  times  it  is  also  desirable  or  necessary  to  look  at 
velocity  vectors.  This  may  be  to  study  the  location  of  separation  and  reattachment 
points  and  to  see  the  size  of  the  separation  zone.  The  development  of  the  bound¬ 
ary  layer  and  shear  layer  can  also  best  be  shown  by  velocity  vector  plots.  The 
present  graphics  code  provides  the  option  of  plotting  velocity  vectors.  To  plot  the 
velocity  vectors,  at  every  grid  point,  a  line  vector  with  or  without  an  arrow  head 
proportional  in  length  to  the  absolute  value  of  the  velocity  at  the  given  grid  point 
is  drawn.  Instead  of  the  velocity  vectors,  an  option  is  provided  to  plot  just  the 
velocity  direction.  In  problems  where  the  magnitude  of  the  velocities  may  change 
drastically,  it  has  been  found  convenient  to  plot  the  velocity  direction. 

Graphics  for  the  Multiple  Grids 

The  multiple  grid  graphics  code  is  based  on  the  single  grid  graphics  code.  In  the 
multiple  patched  grid  solution  procedure  grids  are  constructed  to  overlap.  When 
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more  than  one  grid  is  used  to  solve  the  flow  problem,  the  major  question  that  arises 
is  what  to  do  in  regions  where  the  grids  overlap.  The  solution  is  obtained  in  all 
of  the  grids  and  so  there  are  regions  in  the  flow  field  where  multiple  solutions  are 
available.  The  accuracy  of  the  solution  on  each  grid  is  influenced  by  grid  fineness  and 
grid  topology,  among  overlapping  grids  the  solution  accuracy  can  be  quite  different. 
Since  the  solutions  in  overlapping  regions  do  not  have  the  same  function  values  or 
accuracy,  any  attempt  to  represent  them  there  will  exhibit  some  non-smoothness. 

Though  it  is  true  that  the  solution  in  the  overlapping  region  is  multivalued,  if 
the  solution  procedure  assures  smoothness  and  continuity  to  the  interface  boundary, 
then  any  one  of  the  solutions  in  the  overlapping  region  could  be  used  to  represent 
the  global  solution.  In  the  solution  procedure,  we  order  the  grids  with  assigned 
indices  in  some  sequence.  The  graphics  plotting  is  done  in  the  reverse  sequence.  In 
general,  the  coarsest  grid  is  the  first  grid,  and  any  finer  grid  interior  to  it  will  have 
a  higher  index  value.  To  plot  the  solution  in  the  global  domain,  the  grid  with  a 
highest  number  (finest  grid)  and  most  accurate  solution  will  be  plotted  first.  Next 
the  solution  in  the  lower  grids  will  be  plotted  sequentially.  In  regions  of  overlap, 
only  the  finer  grid  with  higher  index  values  is  used.  It  is  of  interest  to  note  that 
the  base  grid  is  subdivided  by  the  higher  level  grids  and  solution  in  the  finer  grid 
regions  except  for  minimal  necessary  overlap  are  available  only  from  the  finer  grid 
solutions  and  not  from  the  coarse  grid. 

In  Figure  28a  we  have  shown  an  example  of  overlapping  grids  for  a  supersonic 
inlet  problem.  Figure  28b  shows  the  pressure  contour  in  the  global  domain.  To 
obtain  this  plot,  the  pressure  contours  in  the  refined  and  shock  aligned  grid  was 
plotted  first  (with  the  single  grid  contour  plot).  Before  the  contours  in  the  base  grid 
were  plotted,  the  previously  plotted  refined  grid  region  was  “blanked  out”  so  that 
no  contour  lines  of  the  coarse  grid  solution  could  subsequently  be  plotted  inside 
the  region.  Blanking  a  curve  bounded  region  is  accomplished  without  special  effort 
using  a  feature  of  the  DISSPLA  graphics  program.  Finally,  the  base  grid  solution 
is  plotted  only  in  the  non-blanked  regions.  For  all  partitioned  grids,  like  the  base 
grid  in  this  case,  we  have  a  special  subroutine,  that  plots  the  contour  curves  on  the 
integrable  cells  of  the  base  grid.  The  part  of  a  contour  line  outside  of  a  blanked 
region  for  a  cell  partially  overlapping  such  a  region  is  plotted  up  to  the  blanking 
boundary. 

The  above  choice  of  solution  representation  does  not  guarantee  smoothness 
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and  continuity  of  contours  along  the  boundaries  of  overlapping  grids.  The  available 
discrete  solution  itself  is  nonsmooth  due  to  spatial  discretization  error  being  different 
in  each  grid.  Without  adding  additional  smoothing,  it  is  not  possible  to  represent 
the  multiple  grid  solution  with  smoothness  and  continuity.  We  have  chosen  not 
to  add  any  smoothing  but  to  represent  the  best  available  solution.  The  technique 
has  a  virtue  of  exhibiting  the  extent  of  truncation  errors  between  grids  and  making 
evident  places  where  more  refinement  may  be  needed.  In  our  numerical  studies,  the 
above  choice  of  contour  plotting  in  multiple  grids  seems  to  represent  the  solution 
very  well. 

3.  Professional  Personnel 

Professional  researchers  who  contributed  to  this  project  are 

Dr.  Charles  K.  Lombard,  Principal  Investigator 

Professor  Joseph  Oliger,  Consultant 

Dr.  Marcel  Vinokur,  Consultant 

Dr.  Ethiraj  Venkatapathy 

Dr.  Jorge  Bardina 

Dr.  N.  Nagaraj 

4.  Interactions 

The  research  described  in  this  report  has  to  this  point  been  partially  presented 
in  the  form  of  a  paper  on  algebraic  grid  generation  by  Vinokur  and  Lombard  (refer¬ 
ence  3),  a  SIAM  meeting  paper  (reference  25)  by  Oliger  and  Lombard  on  boundary 
procedures  for  bidiagonal  alternating  sweep  schemes  and  a  Computational  Fluid 
Dynamics  Seminar  at  NAS  A- Ames  Research  Center  by  Lombard  on  the  Univer¬ 
sal  Single  Level  Implicit  Algorithm.  The  latter  invoked  tremendous  interest  and 
discussion. 

The  research  as  spawned  three  meeting  papers  on  the  single  level  relaxation 
algorithm  -  references  22,  23  and  24.  Four  papers  have  been  written  on  interior 
patch  boundary  treatment,  data  structure  and  applications  of  patched  mesh  sys¬ 
tems.  There  are  references  27,  28,  30  of  the  present  report  and 

Lombard,  C.K.  and  Venkatapathy,  Ethiraj:  “Implicit  Boundary  Treatment  for 
Joined  and  Disjoint  Patched  Mesh  Systems,”  Workshop  on  High  Reynolds  Flows, 
Nobeyama,  Japan,  September  1985. 
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The  research  is  about  to  give  rise  to  a  paper  on  simulation  of  multi  rocket 
engine  exhaust  flow,  a  problem  of  interest  to  and  partially  supported  by  AFRPL. 

5.  New  Discoveries 

The  Universal  Single  Level  Algorithm  for  the  compressible  Euler  or  Navier- 
Stokes  equations  is  a  new  discovery  in  numerical  methods  that  promises  to  result 
in  substantial  efficiencies  in  data  storage,  programming,  machine  time  and  human 
productivity. 
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Figure  1.  Portion  of  an  al gebra ical ly  generated 
computational  mesh  for  a  flow  domain  containing 
external  and  an  internal  corner. 


gure  Representative  2-D  mesh  generated  with 
mplified  interactive  algebraic  procedure. 
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Figure  3.  Representative  3-D  mesh  generated  with  simplified 
interactive  algebraic  procedure. 


Figure  4.  Shubin's  diverging  nozzle  supersonic 
flov/  solution  developed  in  one  forward  sweep  from 
supersonic  initial  data. 
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in's  diverging  nozzle  flow  solution  developed  in  alternating  forward 
backward  sweeps,  one  each  per  global  iteration  step  for  five  steps, 
re  -  computed  results  ,  line  -  exact  solution 
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Figure  5.  continued 
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Figure  6.  Blottner's  converging-diverging  supercritical  nozzle  flow  solution 
developed  in  alternating  direction  sweeps  for  ten  global  iteration 
steps.  square  -  computed  results  ,  line  -  exact  solution 
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based  on  the  exact  steady  solution  for  the  the  full  block  tridiagonal  CSCM  Solver  and  the 

subcritical  nozzle  flow  solved  with  forward  CSCM-S  method.  The  solutions  have  reached  an  RMS 

sweeps  only.  error  less  than  1  x  10-5. 
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Figure  ] 4 .  Schematics  of  the  supersonic  inlet  flow  problem 
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Figure  17.  Wall  pressure  comparison.  Line  -  exact  inviscid 
solution,  square  -  first  order  inviscid  solution, 
dot  -  first  order  viscous  solution. 
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Figure  24.  Test  on  three  patches  of  implicit 
stability  and  rate  of  convergence  of  (ci rcl es ) 
computed  boundary  point  operator  differencing 
to  frozen  data  in  adjacent  patches;  (solid 
line)  effectively  continuous  grid  method. 


Figure  27.  Pressure  contours 
computed  on  101  x  101  point 
uni  form  mesh. 


27a .  First  order . 


27b .  Second  order . 
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Figure  28a.  Shock  aligned  patched  grids 
for  the  inlet  problem. 


Figure  28b.  Pressure  contours  from  first 
order  solution  on  patched  grid. 


Figure  28d .  Sketch  of  patched  adaptive 
mesh  topology  concluded  from  present 
results  to  be  effective  for  capture  of 
flow  structure  of  the  inlet  problem. 


Figure  28c.  Pressure  contours  from  sec¬ 
ond  order  solution  on  patched  grid. 


